Path Integral Monte Carlo study of phonons in the bcc phase of 4 He 
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Using Path Integral Monte Carlo and the Maximum Entropy method, we calculate the dynamic 
structure factor of solid 4 He in the bcc phase at a finite temperature of T = 1.6 K and a molar 
volume of 21 cm 3 . Both the single-phonon contribution to the dynamic structure factor and the total 
dynamic structure factor are evaluated. From the dynamic structure factor, we obtain the phonon 
dispersion relations along the main crystalline directions, [001], [Oil] and [111]. We calculate both 
the longitudinal and transverse phonon branches. For the latter, no previous simulations exist. We 
discuss the differences between dispersion relations resulting from the single-phonon part vs. the 
total dynamic structure factor. In addition, we evaluate the formation energy of a vacancy. 

PACS numbers: 67.80.-s, 05.10.Ln, 63.20.Dj 



I. INTRODUCTION 

Solid helium, the well-known example of a quantum 
solid, continues to be a subject of interest to theorists 
and experimentalists alike. It is characterized by large 
zero-point motion and significant short-range correlation 
of its atoms. These effects make the theoretical descrip- 
tion of the solid very difficult. The self - consistent 
phonon (SCP) method, 1,2 which has been developed over 
the years to treat this problem, takes into account the 
high anharmonicity and short-range correlations in or- 
der to calculate the dynamical properties of solid helium. 
The predictions of the SCP agree well with experiment 
in the fee and hep solid phases. In the low-density bcc 
phase, the agreement between the theory and experiment 
is less satisfactory, in particular regarding the transverse 
phonons along the [110] direction. 3 The SCP theory is 
a variational perturbative theory, and is implemented at 
zero-temperature. 1 As a complementary approach to the 
SCP, numerical simulations have been performed over 
the years. Boninsegni and Ceperley 4 used Path Integral 
Monte Carlo (PIMC) to calculate the phonon spectrum 
of liquid 4 He at finite temperature. Galli and Reatto 5 
used the Shadow Wave Function approach to obtain the 
spectra of longitudinal phonons in hep and bcc 4 He at 
zero temperature. 

Recently, interest in the properties of quantum solids 
has been revived, following reports indicating the pos- 
sible existence of a "supersolid" in the hep phase 6 . In 
addition, an optic-like excitation branch was recently 
discovered in the bcc phase by Markovich et al, 7,8 (in 
a mono-atomic cubic solid, one should observe only 3 
acoustic phonon branches). These results indicate that 
the physics of solid helium is not yet entirely understood. 
Gov et. al. 3 proposed that the new excitation branch is a 
result of the coupling of transverse phonons to additional 
degrees of freedom, unique to a quantum solid. 

In order to reexamine some of these issues using an al- 
ternative approach, we decided to study the excitations 
in bcc solid helium 4 He, by performing Quantum Monte 
Carlo (QMC) simulations at a finite temperature. We 
use Path Integral Monte Carlo (PIMC) 10-12 , which is a 



non-perturbative numerical method, that allows, in prin- 
ciple, simulations of quantum systems without any as- 
sumptions beyond the Schrodinger equation. The two 
body interatomic He-He potential 13 is the only input for 
the PIMC simulations. In our study the Universal Path 
Integral (UPI) code of Ceperley 10 was adapted to calcu- 
late the phonon branches at finite temperature. 

The novel features of our study include the calculation 
of the transverse phonon branches of bcc 4 He at a finite 
temperature of 1.6 K where this phase is stable. Trans- 
verse phonons are of particular interest due to their pos- 
sible relation with the new optic-like excitation 3 . For lon- 
gitudinal phonons, we observe a difference between dis- 
persion relations resulting from the single-phonon part of 
dynamic structure factor and the total structure factor. 
This difference becomes significant at large wavevectors. 
Finally, we repeated our calculations of the longitudinal 
phonon spectra in the presence of point defects and eval- 
uate the formation energy of a vacancy at a constant 
density and at a constant volume. We describe details 
of our simulations in Sec. II. In Sec. Ill we present the 
results of our calculations, and summarize them in Sec. 
IV. 



II. METHOD 
A. Theory 

The PIMC method used in our simulations is based on 
the formulation of quantum mechanics in terms of path 
integrals. It has been described in detail by Ceperley 10 . 
The method involves mapping of the quantum system of 
particles onto a classical model of interacting "ring poly- 
mers" , whose elements, "beads" or "time-slices" , are con- 
nected by " springs" . The method provides a direct sta- 
tistical evaluation of quantum canonical averages. In ad- 
dition to static properties of the system, dynamical prop- 
erties can be also extracted from PIMC simulations. 10 

The object of this study is the phonon spectrum, which 
can be extracted from the dynamic structure factor, 
S(Q,u>). We would like to express S(Q,ui) in terms of 



2 



phonon operators 1 . The definition of S(Q, u) in terms of 
density fluctuations is 



Fq,\(t) to the complex plane 4 t 
time, we obtain 



it. Using imaginary- 



die™ < p Q (t)p_ Q (0) >, (1) 



where HQ and ftio are the momentum and energy ( we 
take h = 1), pQ is the Fourier transform of the density of 
the solid, and n is the number density. <S(Q, w) is usually 
expressed in terms of phonons, by writing S(Q,lu) as a 
sum of terms involving the excitation of a single phonon, 
Si(Q,u), a pair of phonons, S*2(Q,w) and higher order 
terms which also include interference between different 
terms. 1,2 In most of our simulations we calculated the 
Si(Q,w) term. Some calculations of S(Q,uj) were also 
performed, and will be discussed below. 

Taking the instantaneous position r(Z,t) of atom I 
as the lattice point R ; plus a displacement u(l,t) — 
r(l, t) — R(Z, t), we rewrite S*(Q, u) in terms of these dis- 
placements. The one-phonon contribution is then given 
by 1 

5 1 (Q,c) = d 2 (Q)^ e ^ R '- R °)([Qua,t)][Qu(/,0)]), 
i 

(2) 

where d(Q) =< cxp(-i(uQ) 2 ) > is the Debye-Waller 
factor. The displacement u(l,t) can be expressed using 
the phonon operators Aq y \(t) 14 ' 15 
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Fq,x(r) = J Si, A (q, w) (e- wr + e^"')) dco (7) 

where -^^(r) is the intermediate scattering function, 
and = 1/kT. 

In our simulations, we sampled the displacement 
u(Z, r) for each "time-slice" t of the l-th atom represented 
by a "ring polymer", and calculated A^\(t) by perform- 
ing spatial Fourier transformation 



Aj,a(t) = ^2 e\u(l, t) exp (iqR ( 



(8) 
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where q is the phonon wave- vector, A is the phonon 
branch index, and e\ are polarization vectors, chosen 
along the directions [001], [Oil] and [111]. Using Aq t \(t), 
the one-phonon term S\(Q,u}) for a specific phonon 
branch is rewritten as 1 



Using (8), Fq is obtained as a quantum canonical 
average of the product of the phonon operator < 
Aq,\(T)A^q_\(0) > in equilibrium. 

In order to calculate 5i ;A (q, uj) from Eq.(7), we need 
to perform an inverse Laplace transformation. Perform- 
ing this inversion is a difficult numerical problem, 10 ' 16 
because of the inherent statistical uncertainty of noisy 
PIMC data. The noise rules out an unambiguous recon- 
struction of the 5i,A(q, w). The best route to circumvent 
this problem is to apply the Maximum Entropy (Max- 
Ent) 16,17 method that makes the Laplace inversion bet- 
ter conditioned. 

The MaxEnt method yields a dynamic scattering func- 
tion, Si t \(q,u>) which satisfies Eq. (7) and at the same 
time maximizes the conditional probability imposed by 
our knowledge of the system. This can be done if some 
properties of <Si ; A(q, u>) are known. For example, the dy- 
namic scattering factor is a non - negative function, and 
has certain asymptotic behavior at small and large u). In 
the MaxEnt method the probability to observe of a given 
dynamic scattering function is given by 
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where AQ, q _G is the delta function, and G is a recip- 
rocal lattice vector. We use Q, which lies inside the 
first Brillouin zone and parallel to one of e\. Therefore, 
Si(Q,w) = <Si(q, u), and Si(q, w) is given by 



(9) 

(4) where P is the probability to observe Si,a(Q, u>) for given 
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and 
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set of sampled < Fq t \(r) >, x is the likelihood, a is a 
parameter and S e „t is the entropy 16 . To simplify our no- 
tation, we use Si (to) below to denote by the one-phonon 
dynamic structure Si.,\(q, uj) for a given q and A, and 
omit the explicit dependence on q and A. Similarly, we 
replace F q> ^(r) by F T . The likelihood x 2 is given by 



X 



where the kernel K TyU1 is defined as 

K TyU> = cxp(-Tw) + exp(— (/3 — uj)t) 



jSi(u))— < F T >) 
(10) 

(11) 



is the intermediate scattering function. 

We cannot directly follow the dynamics of helium 
atoms in real time using the Quantum Monte Carlo 
method. However, we can extract information about 
the dynamics by means of the analytical continuation of 



The covariance matrix, C T y, describes the correlation 
between the different time slices r for a given atom 
("ring" polymer). This matrix is defined as 



C' T T > =< F T F T > > - < F T >< F T > >, 



(12) 
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where < F T > is obtained as an average over all atoms at 
a given time slice r. Because F T is periodic as one goes 
around the polymer 4 , the summation on r is done for 
t = 1, M/2, where M is the total number of time-slices 
in a "polymer ring" . 

The entropy term S ent is added to x 2 in order to make 
the reconstruction procedure better conditioned 10,16 . We 
remark here that in some QMC simulations only the diag- 
onal elements of C TjT ' are taken into account 16 , but here 
we use all the elements, because the < F T > at different 
r are correlated with each other. 

Although x" 2 measures how closely any form of S\(u>) 
approximates the solution of Eq. (7), one cannot deter- 
mine Si(u>) reliably from PIMC using y 2 alone. 10,16 To 
make this determination, one needs to add the entropy 
term S en t to x 2 in (9) to make the reconstruction proce- 
dure better conditioned. The entropy term is given by 




co(K) 



Sent(w) = - [ da; ( Si(oj) In ^1^1 + Si(w) - m(oj) J 
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(13) 

where m{u>) includes our prior knowledge about the prop- 
erties of Si(lo), examples of which were given above. 
The simplest choice of S%(u>) is the fiat model, in which 
m(u;)=const for a selected range of frequencies and zero 
otherwise. We took a cutoff frequency corresponding to 
an energy of 100K. This flat model is used as an input 
for most of our simulations. In addition, we used a "self- 
consistent" model where the output of a previous Max- 
Ent reconstruction is used as input for the next MaxEnt 
reconstruction in an iterative fashion 4,10 The flat model 
is taken for the initial iteration. Finally, we also tried 
m(uj) with Gaussian and Lorentzian shape, with peaks 
given by the SCP theory and experiment. As explained 
below, the outcome is not very sensitive to the choice of 
m(u) provided the PIMC data are of good quality. 

Finally, we discuss the parameter a in (9). The mag- 
nitude of this parameter controls the relative weight of 
the PIMC data vs. the entropy term in the determina- 
tion of Si(lui). There are different strategies to obtain a. 
In our simulations we used both the "classical" MaxEnt 
method 16 and random walk sampling. 4 The "classical" 
MaxEnt method picks the best value of a, while random 
walk sampling calculates a distribution of a, 7r(a). Next, 
for each value of a we calculate S%(u>). The final Si(w) 
is obtained as a weighted average over a. We found that 
when the collected PIMC data is of good quality, the dis- 
tribution 7r(a) becomes sharply peaked and symmetric, 
and the phonon spectra obtained by means of "classical" 
MaxEnt and random walk are almost the same. Good 
quality data are characterized by an absence of correla- 
tion between sequential PIMC steps and by small statis- 
tical errors. 



FIG. 1: (Color online) Longitudinal component of the dy- 
namic structure factor, S(q,uj), for q = 0.83 r.l.u. along the 
[100] direction: single phonon contribution (circles, red on- 
line), and total structure factor (squares, black online). The 
lines are a guide to the eye. 



B. Our implementation 

The MaxEnt method assumes that the distribution of 
the sampled F T is Gaussian. 16 We re-block 16 the sam- 
pled values of F T in order to reduce the correlations and 
to make the distribution as close as possible to a Gaus- 
sian, with zero third (skewness) and fourth (kurtosis) mo- 
ments. 

The criterion determining the minimum number of 
sampled data points comes from the properties of the 
covariance matrix C T)T <. If there are not enough blocks 
of data the covariance matrix becomes pathological. 16 
Therefore, the number of blocks must be larger than the 
number of time slices in a "ring" polymer. In our sim- 
ulations, each atom is represented by a "ring" polymer 
with 64 time slices. We collected at least 300 blocks in 
each simulation run. We found that at least 10000 data 
points were required in order to obtain the 300 blocks. 
Each simulation run took about two weeks of 12 Pentium 
III PCs running in parallel. 

Statistical errors were estimated by running the PIMC 
simulations 10 times, with different initial conditions in 
each case. After each run, S(q, to) was extracted using 
the MaxEnt method. The phonon energy for a given 
q was then calculated by averaging the positions of the 
peak of S(q,u) over the set of the simulation runs. The 
error bars of each point shown in the figures below rep- 
resent the standard deviation. 

In the simulations we used samples containing between 
128 and 432 atoms. This allowed us to calculate 5(q, ui) 



for values of q between 0.17 and 1 in relative lattice units 
(r.l.u.= 2n/a, where a is the lattice parameter). The 
number density was set to p = 0.02854 (1/A 3 ) and the 
temperature was T = 1.6 K. A perfect bcc lattice was 
chosen as the initial configuration. The effects of Bose 
statistics are not taken into account in our simulation, 
which is a reasonable approximation for the solid phase. 
A typical example of the calculated dynamic structure 
factor is shown in Fig.l. The figure shows both the sin- 
gle phonon contribution Si(q, id) and the total S(q, id) 
for a longitudinal phonon along the [001] direction. To 
illustrate the difference between Si(q, id) and S(q, id), we 
chose to show the results for q close to the boundary of 
the Brillouin zone. This difference is discussed below. 



III. RESULTS 
A. Phonon spectra 

The calculated longitudinal and transverse phonon 
spectra of solid 4 He in the bcc phase along the main 
crystal directions ([001], [111] and [011]) are shown in 
Figs. 2-7. We compare our results with the experimen- 
tal data measured by inelastic neutron scattering from 
bcc 4 He with a molar volume of 21.1 cm 3 at T = 1.6 K., 
by Osgood et al, 18-20 and by Markovitch et al. 8 




0.2 0.4 0.6 0.8 1 
q (r.l.u) 

FIG. 2: (Color online) Calculated dispersion relation of the 
L[001] phonon branch (squares, red online) using Si(q,(d). 
Experimental data are from 18,20 (triangles up) and from 8 (tri- 
angles down) . The error bars represent statistical uncertainty. 

As expected, the agreement between our simulations of 
Si (q, uj) and experiment is very good at small q, where 
one-phonon excitation is the most significant contribu- 
tion to S^q, uj). As q increases, higher order processes be- 



4 




Q| , 1 , 1 , 1 , 1 , 

0.2 0.4 0.6 0.8 1 
q (r.l.u.) 

FIG. 3: (Color online) Calculated dispersion relation of the 
T[001] phonon branch (squares, red online) using Si(q,ai). 
Experimental data are from 18 ' 20 (triangles up) and from 8 (tri- 
angles down) . The error bars represent statistical uncertainty. 



T 




q (r.l.u.) 

FIG. 4: (Color online) Calculated dispersion relation of the 
L[011] phonon branch (squares, red online) using 5i(q, id). 
Experimental data are from 18,20 (triangles up) and from 8 (tri- 
angles down) . The error bars represent statistical uncertainty. 



come significant, and the calculated values deviate from 
the experimental data, especially along [001] and [111]. 
In the case of longitudinal phonons, it is possible to cal- 
culate their energies using the total <S(q, id) obtained di- 
rectly from Eq. (1) instead of just the single phonon 
contribution. The dispersion relations calculated with 
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Q| , 1 , 1 1 1 , 1 1 |_| 

0.1 0.2 0.3 0.4 0.5 
q (r.l.u.) 

FIG. 5: (Color online) Calculated dispersion relations of 
transverse phonon branches along [011] using Si(q,a>). Calcu- 
lated values are shown for the Ti branch (squares, red online) 
and T2 branch (diamonds, blue online). Experimental data 
are from 18,20 ( Ti-triangles up, T2-triangles left) and from 8 
(Ti-triangles down, T2-triangles right). The error bars repre- 
sent statistical uncertainty. 



S{<\,u>) are shown in Figs. 8 and 9. It is evident that 
using the total scattering function improves the agree- 
ment with experiment at large q, especially for the [111] 
direction. 

We point out that the calculated phonon branches 
Ti(110),T 2 (110) and L(110) show good fit to the exper- 
imental data. Our results were obtained with the two 
body potential, which takes the He atoms as point parti- 
cles. Gov et al. 3 suggested that one needs to go beyond 
this approximation to obtain the Ti(llO) phonon branch 
in good agreement with experiment. Gov's approach also 
predicts the new excitation branch observed recently 7 . 
Although the calculated Ti(110) branch is in agreement 
with experiment without any additional assumptions, we 
were not able to see the new excitation in our simula- 
tions. Experimentally, this excitation is about an order 
of magnitude less intense than a phonon. It is best ob- 
served in scattering experiments done with very small 
q < 0.1 r.l.u 8 . Both these factors make it very difficult 
to search for this excitation in simulations. Whether it 
can be found in this approach remains an open question. 

In addition to experimental results, our simulations 
can also be compared with those of Galli and Reatto 5 , 
who used the Shadow Wave Function (SWF) approach to 
calculate the longitudinal phonon branches of bcc 4 He. 
As shown in Fig. 10, the overall agreement between these 
PIMC simulations and SWF results is good. 
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0.2 0.4 0.6 0.8 1 
q (r.l.u.) 

FIG. 6: (Color online) Calculated dispersion relation of the 
L[lll] phonon branch (squares, red online) using Si(q,u). 
Experimental data are from 18,20 (triangles up) and from 8 (tri- 
angles down). The error bars represent statistical uncertainty. 




0' 1 1 1 1 1 1 ' 1 ' 

0.2 0.4 0.6 0.8 

q (r.l.u.) 

FIG. 7: (Color online) Calculated dispersion relation of the 
T[lll] phonon branch (squares, red online) using Si(q, k>). 
Experimental data are from 18,20 (triangles up) and from 8 (tri- 
angles down) . The error bars represent statistical uncertainty. 



B. Vacancies 

Recent experimental work 6 revived the interest in 
point defects, such as vacancies. It is therefore inter- 
esting to examine the influence of vacancies on the prop- 
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FIG. 8: (Color online) Calculated dispersion relation of 
the L[001] phonon branch (squares, red online) using 
S'(q, a;). Experimental data are from 18,20 (triangles up) and 
from 8 (triangles down). The error bars represent statistical 
uncertainty. 



FIG. 9: (Color online) Calculated dispersion relation of 
the L[lll] phonon branch (squares, red online) using 
S(q,uj). Experimental data are from 18,20 (triangles up) and 
from 8 (triangles down). The error bars represent statistical 
uncertainty. 



erties of the solid. We repeated our calculation of the 
phonon branches in the presence of 0.23% vacancies (1 
atom of 432). Within the statistical error bars, we found 
no difference between the phonon energies with or with- 
out vacancies. Galli and Reatto 5 found that vacancies 
lower the energies of the phonons close to the boundary 
of the Brillouin zone. However, in their simulation they 
used a concentration of vacancies of 0.8%, so that the 
cumulative effect may be larger. We also calculated the 
vacancy formation energy, AE V , according to Pederiva et 
al. 21 

AE v = (E(N-l,p)-E(N,p))(N-l), (14) 

where E(N, p) is the total energy of N atoms. The en- 
ergy E(N,p)) was calculated for a perfect crystal, while 
E(N — l,p) was calculated after removing one atom. 
The density of two systems was kept the same by ad- 
justing the lattice parameter. Values of AE V calculated 
using the PIMC, Shadow Wave Function (SWF) 9 and 
Shadow Path Integral Ground State (SPIGS) 5 methods 
are summarized in Table F In addition, we calculated 
AE V at constant volume, which is the condition usu- 
ally realized in experiments rather than constant density. 
We obtained AE V = 5.7 ± 0.7 K. The lower value arises 
since the repulsive part of the potential is weaker in a 
sample having lower density. There is no generally ac- 
cepted experimental value 22 of AE V . According to NMR 
studies 23- the energy of vacancy formation in the bec 
phase is AE V = 6.5 ± 0.2(K), while X-ray studies 25 sug- 
gest that AE V = 9±1(K). We comment here that the cal- 
culated values of AE V are significantly smaller than 14K, 



TABLE I: Calculated energy of formation of a vacancy, AE V , 
for bec solid 4 He. N is the number of atoms used in each of 
the simulations. 



source 


method 


density {l/A 6 ) 


N 


AE V (K) 


this work 


PIMC 


0.02854 


128 


10.57 ±0.38 


this work 


PIMC 


0.02854 


250 


9.96 ± 0.89 


ref. y 


SWF 


0.02854 


128 


8.08 ±2.76 


ref. y 


SWF 


0.02854 


250 


6.69 ± 3.86 


ref. b 


SWF 


0.02898 


128 


8.9 ±0.3 


ref. 5 


SPIGS 


0.02898 


128 


8.0 ±1.3 



the energy of the new excitation observed by Markovitch 
et al. 7 ' 8 . Hence, this new excitation does not seem to be 
a simple vacancy. 



IV. CONCLUSIONS 

We calculated the dynamic structure factor for solid 
helium in the bec phase using PIMC simulations and 
the MaxEnt method. PIMC was used to calculate the 
intermediate scattering function in the imaginary time 
from which the dynamic structure factor was inferred 
with the MaxEnt method. We extracted the longitudi- 
nal and transverse phonon branches from the one-phonon 
dynamic structure factor. To the best of our knowledge 
this is the first simulation undertaken for the transverse 
branches. At small q, where the one-phonon excitation 
is the most significant contribution to the dynamic struc- 
ture factor, the agreement between our simulations and 
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0.5 
q (r.l.u.) 



FIG. 10: (Color online) A comparison ol dispersion relations 
of the L[001] phonon branch obtained in the present work us- 
ing S(q, lu) (squares, red online), with the same relation cal- 
culated by means of the Shadow Wave Function technique 5 
(circles). The error bars represent statistical uncertainty. 



experiment is very good. At large q, multi-phonon scat- 
tering and interference effects 1 becomes important. Con- 
sequently, the position of the peak of in the S\ (q, u) does 
not correspond to the position of the peak in the 5'(q, lu), 
and the phonon energies calculated from Si (q, uj) are too 
low. If S(q,iv) is used instead of Si(q, to), the agreement 
with experiment is significantly improved. We repeated 
the simulations in the presence of 0.23% of vacancies, and 
found no significant differences in the phonon dispersion 
relations. We also calculated the formation energy of a 
vacancy both at constant density and at a constant vol- 
ume. 
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